Theoretical Support for the Hydrodynamic Mechanism of Pulsar Kicks 



J. Nordhaus/'Q T. D. Brandt, 1 A. Burrows, 1 E. Livne, 2 and C. D. Ott 3 

'Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, U.S.A. 
2 Racah Institute of Physics, Hebrew University, Jerusalem, Israel 
'Theoretical Astrophysics, Mail Code 350-17, California Institute of Technology, Pasadena, CA 91125 U.S.A. 

(Dated: November 8, 2010) 

The collapse of a massive star's core, followed by a neutrino-driven, asymmetric supernova explo- 
sion, can naturally lead to pulsar recoils and neutron star kicks. Here, we present a two-dimensional, 
radiation-hydrodynamic simulation in which core collapse leads to significant acceleration of a fully- 
formed, nascent neutron star (NS) via an induced, neutrino-driven explosion. During the explosion, 
a ~10% anisotropy in the low-mass, high-velocity ejecta lead to recoil of the high-mass neutron star. 
At the end of our simulation, the NS has achieved a velocity of ~150 kms -1 and is accelerating 
at ~350 kms -2 , but has yet to reach the ballistic regime. The recoil is due almost entirely to hy- 
drodynamical processes, with anisotropic neutrino emission contributing less than 2% to the overall 
kick magnitude. Since the observed distribution of neutron star kick velocities peaks at ~300-400 
kms -1 , recoil due to anisotropic core-collapse supernovae provides a natural, non-exotic mechanism 
with which to obtain neutron star kicks. 

PACS numbers: 97.60.Bw, 97.60.Gb, 97.60.Jd, 95.30.Jx, 95.30.Lz 



I. INTRODUCTION 

The velocity distribution of young pulsars bears lit- 
tle resemblance to that of their massive star progeni- 
tors [I]. Typical birth velocities range from ^200-500 
kms -1 , with some reaching upwards of ^1000 kms -1 
[2J. While the observed pulsar velocities may hint at a 
two-component distribution (possibly implying two pop- 
ulations) [SHS], recent work supports a single, Maxwellian 
distribution [gHTUj. 

Various mechanisms for the origin of neutron star kicks 
and pulsar recoil and their connections with pulsar spins 
have been proposed Misaligned jet/counter-jets dur- 
ing the supernova explosion might produce sufficient ac- 
celeration if they are launched near the proto-neutron 
star (PNS) [T2J 03]. However, such jets are generated 
only in fast rotators and may not be generic [TlrfT^] . An- 
other possibility is anisotropic neutrino emission from the 
cooling proto-neutron star. If strong magnetic fields are 
present, neutrino- matter interactions can generate dipole 
asymmetries of ~1%, leading to recoil on the order of 
a few hundred kms -1 P"?rt20| . These scenarios require 
magnetar field strengths (i.e. 10 14 — 10 15 G) and/or ex- 
otic neutrino physics [2"Trf2~i] and may not produce sub- 
stantial kicks in typical core-collapse supernovae. 

If neutron star kicks are a generic feature of core col- 
lapse, then the most natural explanation is recoil due 
to an asymmetric supernova explosion [2"5rf2"8] . During 
axisymmetric core collapse, the stalled bounce shock is 
unstable to neutrino-driven convection and low-order t- 
modes. Significant asymmetry at the onset of neutrino- 



driven shock revival should naturally lead to an asym- 
metric explosion and the hydrodynamic recoil of the PNS 
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Observations of large-scale asymmetries in young su- 
pernova remnants lend qualitative support to the hy- 
drodynamic mechanism |32j . Unfortunately, multi- 
dimensional, radiation-hydrodynamic simulations of re- 
coil are computationally challenging. A proper study 
requires simulating the full physics of collapse, the for- 
mation of the PNS, the development of instabilities dur- 
ing the post-bounce phase, the evolution of the asym- 
metric explosion, the off- axis movement of PNS, and the 
full decoupling of the ejecta from the PNS. Because the 
expanding post-shock material interacts with the PNS 
through both pressure and gravity, this requires following 
the shock out to large distances (hundreds of thousands 
of kilometers) and late times (several seconds) . Compli- 
cating matters is that during this evolution, one must 
continue to resolve the movement of the PNS and the 
surrounding highly nonlinear flow. 

| S check et al. 2006 present a practical approach to this 
problem [251 [55] • By excising the PNS and replacing 
it with a rigid, contracting boundary, they avoid severe 
Courant timestep restrictions. They also greatly simplify 
their radiation transport, enforcing a constant luminos- 
ity at their inner boundary, and begin their calculations 



20 ms after bounce. These approximations allow Scheck 
|et al.| to follow the evolution of the shock to large dis- 
tances and late times, and to perform a detailed param- 
eter study. Unfortunately, this approach requires them 
to infer a kick through a rigid, impenetrable boundary. 
Their results should therefore be checked by more realis- 
tic (though costly) simulations. 



As a complement to the work of Scheck et al. 



we 



present a two-dimensional (2D) simulation of the collapse 
of a 15-Mq progenitor core. By employing a pseudo- 
Cartesian mesh at the center of our domain, we naturally 



2 



capture the neutron star's formation and any subsequent 
off-center acceleration. During our simulation, the proto- 
neutron star forms, after which it recoils due to a delayed, 
neutrino-driven, anisotropic explosion. The explosion is 
artificially induced by adding additional neutrino lumi- 
nosity to the calculation. At the end of our simulation, 
the NS has achieved a velocity of ~150 km s _1 and is 
still accelerating at ~350 kms~ 2 . The recoil is primarily 
hydrodynamic in nature, with anisotropic neutrino emis- 
sion contributing less than 2% of the overall kick magni- 
tude. Most notably, we obtain a significant kick without 
invoking strong magnetic fields, exotic neutrino physics, 
or misaligned jets. Our results are consistent with the 
previous Scheck et al. studies [26l [28]. Taken together, 
these simulations provide compelling numerical support 
for the hydrodynamic mechanism of neutron star kicks. 
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II. 



NUMERICAL SETUP AND METHODS 



Our 2D, axisymmetric calculations are performed with 
the multi-group, arbitrary Lagrangian-Eulerian (ALE), 
radiation-hydrodynamics code VULCAN/2D [33] ■ We 
perform 2D radiation transport using the multi-group 
flux-limited diffusion approximation [34] . We simulate 
the collapse of the inner 5000 km of a non-rotating, 
15- Mq, solar- metallicity, red-supergiant progenitor |35j . 
Exterior to 20 km, our computational domain is a 
spherical-polar mesh which transitions to a pseudo- 
Cartesian grid in the center. Such a grid avoids severe 
timestep restrictions due to the convergence of angular 
zones and frees the PNS to move in response to radia- 
tion or hydrodynamic forces. Our mesh covers the full 
180°, 2D domain with 120 angular zones and 330 radial 
zones (logarithmically spaced exterior to the inner Carte- 
sian region). We employ the finite-temperature nuclear 
equation of state of Shen et al. [361 EI] an d include self- 



gravity with a grid-based solution of the Poisson equation 
[55] . To ensure that we optimally resolve the high-density 
core, we allow our grid to track the PNS. Our remapping 
scheme determines the center of mass of the inner core 
(i.e. densities above 10 12 g cm -3 ) after each timestep and 
shifts the mesh to keep the core centered while ensuring 
momentum conservation. 

Despite decades of intense theoretical effort, the suc- 
cess of the delayed-neutrino mechanism [39 41 in driv- 
ing core-collapse supernova explosions has still not been 
demonstrated [3T1 l4"2"r[5"l"] . However, recent calculations 
have shown that this mechanism's capacity to power 
explosions increases with dimension |52[ 153) . Ambi- 
tious three-dimensional calculations with accurate neu- 
trino transport may yet validate the delayed-neutrino 
mechanism. 

Because previous core-collapse studies with VUL- 
CAN/20 did not produce neutrino-driven supernovae 
[l4l [27] I54H56] , we induce explosions by supplementing 
the radiation transport with additional electron and anti- 
electron neutrino luminosity (L Ve = L Ve — 2 x 10 52 
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FIG. 1: The recoil of the proto-neutron star due to an asym- 
metric core-collapse supernova explosion. The large-scale ex- 
plosion is primarily in the +Z direction (top) while the PNS 
is kicked in the — Z direction (bottom). In the bottom panel, 
the white line denotes Z = 0, while the purple and black 
curves represent the isodensity surfaces where p — 10 12 and 
10 14 g cm" 3 , 
black. 



respectively. Velocity vectors are overlaid in 



ergs -1 ) as described in [52] [53]. This represents an en- 
hancement in the v e and V e luminosities of ^50%. The 
core collapses to nuclear densities, launching a bounce 
shock which stalls and is subsequently revived mainly by 
charged-current neutrino absorption after a delay of ap- 
proximately 135 milliseconds. 
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III. RECOIL FROM ASYMMETRIC 
CORE-COLLAPSE EXPLOSIONS 

At the onset of explosion, the hydrodynamic flow be- 
hind the shock is turbulent and the shock itself is de- 
formed by the development of low-mode instabilities 
[261 EH S3 Ml [58]. The PNS recoils due to the blast's 
anisotropic propagation through the stellar envelope. We 
follow the explosion and the acceleration of the PNS un- 
til 470 ms after bounce, at which point the shock front 
reaches the boundary of our computational domain (5000 
km). Figure [I] shows the global explosion geometry and 
the position of the PNS at the end of our calculation. 
The top panel is an entropy map of our computational 
domain with velocity vectors overlaid and the shock po- 
sition outlined in white. The bottom panel shows the 
electron fraction Y e over the inner ~70 km. The white 
line is the Z = axis, while the pink and black curves 
represent the 10 12 g cm -3 and 10 14 g cm -3 isodensity 
contours, respectively. Note that the asymmetry of the 
explosion in the +Z-direction leads to a PNS recoil in 
the — ^-direction. While axisymmetry restricts our core 
to motion along the Z-axis, three-dimensional computa- 
tions would impose no such constraint and could produce 
a recoil in any direction for initially non-rotating progen- 
itors. Note that the presence of rotation may lead to a 
preferred explosion direction and, hence, kick direction. 
The differences between kicks from non-rotating and ro- 
tating progenitor models should be investigated in 3D. 

While VULCAN/2D automatically and self- 
consistently computes the acceleration of the core, 
it does not output the individual forces governing the 
motion of the PNS. We therefore post-process our results 
by computing the hydrodynamic acceleration a c of the 
core due to anisotropic gravitational forces, pressure 
forces, and momentum flux. The Eulerian equations of 
hydrodynamics give 



Gr , 1 

-^r-am 



r>r c 



PdS + (p pv r vdS 

'c ^ T—r c . 

where p is the density, M c and v c are the mass and mean 
velocity of the inner region (where p > 10 12 gem -1 ), P 
is the gas pressure, v is the fluid velocity, v r is the radial 
component of the velocity, and r c is a fiducial spherical 
radius. The code self-consistently yields the recoil speed 
of the PNS (approximately bounded by the purple curve 
in Fig. [I]), but we can use Eq. [I] to determine the various 
contributions to its acceleration and consequent motion. 

The first term in Eq. [T] represents the acceleration due 
to the gravitational field exterior to r c , assuming a spher- 
ically symmetric distribution of matter interior to this 
radius. The second term is due to anisotropic gas pres- 
sure, while the third term represents the contribution due 
to momentum flux. In a spherically symmetric explo- 
sion, each term would vanish individually. These three 
terms include all hydrodynamic forces, but do not in- 
clude asymmetries in the radiation pressure. In our sim- 
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FIG. 2: The core velocity as a function of time after bounce. 
The solid-red curve in both figures shows the core velocity, 
in the —Z direction, as a function of time after bounce in 
our simulation. Though the inferred core velocities calcu- 
lated at r c = 200 km (top figure) and r c = 500 km (bottom 
figure) accurately reproduce the actual core velocity at late 
times, this figure demonstrates that one must exercise caution 
when interpreting the relative contribution of each compo- 
nent. Anisotropic neutrino flux contributes very little (< 2%) 
of the total kick at all radii. 



ulations, exterior to the radius at which the flux limiter 
transitions to free-streaming, anisotropic neutrino mo- 
mentum contributes ~2% of the total kick (see Fig. [2]). 

In general, the relative contributions of the various 
terms in Eq. [T] will depend sensitively on the radiation- 
hydrodynamics and explosion dynamics. For instance, a 
spherically-symmetric distribution of ejected mass pos- 
sessing asymmetric ejection velocities will lead to gravity 
and momentum terms of the same sign. In particular, 
since the PNS recoils towards the lower-velocity ejecta, 
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the gravitational acceleration is in the same direction as 
the kick. This gravitational "tug-boat" effect enhances 
the recoil. Isotropic ejection velocities with anisotropic 
mass loss results in the gravity component partially can- 
celing the momentum contribution. 

We present the PNS kick velocity (as computed by 
VULCAN/2D) as a solid red line in both the top and 
bottom panels of Fig. [2j Using Eq. [T] we show the in- 
ferred kick velocity (dashed-blue curve) and its compo- 
nents at 200 km (top panel) and 500 km (bottom panel). 
These curves represent the mean velocities of matter in- 
terior to 200 km and 500 km. As the core evolves, matter 
interior to 500 km becomes more centrally concentrated 
and its average velocity approaches that of the innermost 
regions (i.e. the monopole approximation gets better and 
better). The agreement between the red line and the blue 
lines therefore improves with time. 

Figure [2] demonstrates that the kick imparted to the 
PNS may be inferred by evaluating Eq. [I] even at large 
radii. However, the relative contributions of the three 
terms in Eq. [T] differ dramatically. At r c = 200 km, the 
late time evolution of our simulation is dominated by the 
gravitational component, while the momentum and pres- 
sure contributions are of opposite sign and comparable in 
magnitude. For r c = 500 km, the pressure and momen- 
tum contributions are approximately equal (in both sign 
and magnitude) and nearly constant between ^200 ms 
and ^470 ms. The secular evolution of the PNS velocity 
at the end of our calculation is governed by the grav- 
itational component. The one component which does 
not depend strongly on radius is the contribution from 
anisotropic neutrino emission, which is small (<2% of the 
kick). 

The interpretation of the kick (though not its value) 
thus depends on the radius at which the terms of Eq.[l]are 
evaluated. At large radii, pressure and gravity vanish and 
an observer will attribute the entire kick to anisotropic 
momentum flux. The story is very different near the 
PNS itself. Because the inner core is nearly in hydro- 
static equilibrium, pressure and gravity are both very 
large and in balance. An observer in this region would 
remark on the near cancellation of the gravitational and 
pressure terms in Eq. [T] For example, in our calcula- 
tions, with a radius r c that moves inward to always en- 
close 1.3 Mq, these two components of the kick cancel 
to one part in 10 2 . Our results demonstrate the limita- 
tions of interpreting the individual components of Eq. [T] 
Since pressure and gravity do work on expanding matter, 
their contributions to the acceleration decrease in mag- 
nitude relative to the contribution due to the anisotropic 
momentum flux. 



A. Extrapolating the Kick 

Figure [2] indicates that our PNS is still accelerating 
at ^350 kms™ 2 when the shock has reached the bound- 
ary of our computational domain. However, the ejecta 



have not yet decoupled from the core to reach the ballis- 
tic regime. The spatial distributions of momentum and 
velocity offer a hint of the core's future evolution, but 
unfortunately do not permit a straightforward extrap- 
olation. Ideally (though at considerable computational 
expense), this would be handled by remapping our re- 
sults onto a larger grid and continuing a full radiation- 
hydrodynamic calculation. However, momentum and ve- 
locity maps, which we show in Fig. |3j offer a useful pic- 
ture of the ejecta at the end of our calculation. 

The top panel of Fig. [3] shows the velocity of matter 
throughout our computational domain in units of the lo- 
cal escape speed, calculated assuming a spherically sym- 
metric distribution of matter. Because the potential is 
dominated by the PNS, this approximation is extremely 
accurate. The map clearly shows that our model has 
not yet reached the ballistic regime, and that the matter 
behind the shock is still accelerating and evolving dy- 
namically. A significant region of matter at Z ~ —1000 
km seems likely to fall back, while a pocket of material 
at Z ~ 2500 km is expanding at nearly twice the local 
escape speed. The infalling region has only ^20% of the 
momentum in the core and, thus, is unlikely to signifi- 
cantly affect our inferred kick. However, the complexity 
of the hydrodynamics makes it impossible to extrapolate 
by assuming, for example, self-similar expansion. 

The lower panel of Fig. [3] shows the projected Z- 
momentum density, pz = irRpvz- The factor ttR, where 
R is the cylindrical radius, is the length of a semicircle 
of revolution. This projects the half-cylinder defined by 
< < 7r in 3D onto the half-plane X > in 2D, 
so that / pz dX dZ gives the correct value for the total 
Z-momentum. This map shows that the high-velocity 
bubbles at Z ~ 2500 km are regions of low density; 
most of the momentum is concentrated behind the shock 
and in the regions behind the highest velocity ejecta at 
Z — 1000 km. At the end of our calculation the PNS 
is still injecting mass and momentum into these regions. 
There appears to be no such injection of momentum into 
the regions at negative Z. If this causes the expansion of 
matter to slow in the —Z direction, it could help maintain 
an asymmetric matter distribution, and thus the gravita- 
tional component of its acceleration, for several seconds. 

The continued acceleration of the PNS will depend 
on the evolution of the asymmetry of shocked material. 
There are a variety of ways to quantify this asymmetry, 
as discussed in [SHIEHj. We choose a = (v z )/{\v\), where 
() denotes a mass-weighted average over the post-shock 
region with r > 100 km (to exclude the PNS itself). This 
is similar to the a presented in |28] . If we assume this 
asymmetry to be constant in time, material on one side 
of the PNS will be a factor of 1 — a as close as material 
on the other side. We may then crudely estimate the 
gravitational acceleration of the core, a c grav as 
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for small a, where r s h is the shock radius and M sb is 
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the total mass of ejecta and shocked envelope material. 
In our calculation, a ~ 0.1 from 300 milliseconds to 470 
milliseconds after bounce. Assuming M s h ~ Af Q and 



a ~ 0.1, then for a. 



c,grav 



to be of order 1 km 



-2 



we 



need to follow the shock out to ~10 5 km. This corre- 
sponds to 5 seconds at a shock velocity of 20,000 kms -1 , 
and represents a challenging computational problem. We 
hope ultimately to address this problem with CASTRO 
[531 [53], a new adaptive mesh refinement radiation- 
hydrodynamics code, which will allow us to follow the 
shock while still resolving the PNS. 



B. Comparison to Previous Work 

Our approach of following the collapse of a massive 
star's core, the formation of a natal PNS, and the subse- 
quent off-axis motion complements previous studies that 
infer kicks on an excised PNS [25]. By omitting 
the inner regions, starting the simulation ^20 ms after 
bounce, and imposing a constant inner neutrino luminos- 
ity, | S check et al.| greatly reduced the problem's compu- 
tational cost. They were thus able to follow the shock 
evolution to large distances (> 10 4 km) and late times 
(>1 s). To approximate a physical neutron star, those 
authors used a contracting inner boundary motivated by 
radiation- hydrodynamic simulations [25] • While attrac- 
tive for calculating long-term evolution, their approach 
requires one to infer a PNS kick through a rigid bound- 
ary of infinite inertial mass. This assumption neglects 
effects resulting from displacement of the PNS relative 
to the surrounding fluid. To compensate, in a subset 
of their simulations, these authors artificially add the 
inferred kick velocity to the gas, mimicking movement 
of the PNS. Our work handles all of these effects self- 
consistently, providing an important check on the various 
approximations made in |26[ 128] . 

Another difference between our work and that of 
|Scheck et al.| is that we implement the momentum equa- 
tion in conservative form using a grid-based solution 
to the Poisson equation. As a result, our model con- 
serves total momentum to better than 1% of the core's 
final value. |S check et al.| solve the Poisson equation us- 
ing a Legendre expansion with a relativistic correction 
[251 [BUI loT] . Recently, Wongwathanarat et al. performed 
a three-dimensional study using the same techniques in 
the IScheck et al.l two-dimensional studies and arrived at 
similar conclusions. 

Given the differences in our complementary tech- 
niques, the agreement between our results and those of 
Scheck et al. is gratifying. Our detailed calculations of 



the first few hundred milliseconds including the core sup- 
port the work of [2(3 I25J > while their extended calcula- 
tions indicate that a final kick magnitude of at least 400- 
500 kms -1 may be likely for our model. Taken together, 
this body of work strongly supports the case that asym- 
metric supernova explosions lead naturally to substantial 
recoil of the PNS. 
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FIG. 3: Top: The ratio of the fluid velocity, v, to escape speed, 
^csc , as a function of position 470 ms after bounce. Bottom: 
Projected Z-momentum density pz as a function of time 470 
ms after bounce. The cylindrical volume element is included, 
so that j pz dX dZ gives the total Z-momentum. We have 
overlaid velocity vectors and a thick black curve representing 
the position of the shock on both panels. 



IV. CONCLUSIONS 

In this work, we have presented the first multi- 
dimensional, multi-neutrino-energy-group, radiation- 
hydrodynamic simulation of a core-collapse supernova 
that results in a formation and acceleration of a nascent 
neutron star. The recoil of the PNS naturally arises from 
the asymmetric nature of the neutrino-driven explosion. 
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At the end of our simulation the PNS has reached a ve- 
locity of ~150 km s _1 , but is still accelerating at ^350 
kms~ 2 . While it is difficult to extrapolate the acceler- 
ation to later times, our PNS would need to maintain 
this value for only a few hundred milliseconds more to 
reach the peak of the observed pulsar velocity distribu- 
tion. This is suggested by Fig. [3] the continued ejection 
of momentum in the + indirection could maintain the 
asymmetric matter distribution and continue to gravita- 
tionally accelerate our PNS. It should also be noted that 
the highest observed kicks (those upwards of 1000 km 
s -1 ) may result from the most asymmetric and energetic 
explosions. 

Hydrodynamic recoil due to neutrino-driven, core- 
collapse supernovae provides a natural mechanism for 
accelerating neutron stars and pulsars without the need 
to appeal to anisotropic neutrino emission or more ex- 
otic scenarios. However, a definitive confirmation of this 
mechanism will require a self-consistent model of core- 
collapse supernova explosions. To avoid constraints im- 
posed by axisymmetry, future work should investigate 
recoil and explosion anisotropies in three dimensions and 
compare the resulting kick velocities with observations. 
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